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Abstract. - In this work, we study the synchronization of coupled phase oscillators on 
the underlying topology of scale-free networks. In particular, we assume that each network's 
component is an oscillator and that each interacts with the others following the Kuramoto 
model. We then study the onset of global phase synchronization and fully characterize the 
system's dynamics. We also found that the resynchronization time of a perturbed node decays 
as a power law of its connectivity, providing a simple analytical explanation to this interesting 
behavior. 



The behavior of an isolated generic dynamical system in the long-term limit could be 
described by stable fixed points, limit cycles or chaotic attractors [1]. We have also learned in 
recent years that when many of such dynamical systems are coupled together, new collective 
phenomena emerge. In this way, the study of regular networks of dynamical systems have shed 
light on a number of natural phenomena ranging from earthquakes to ecosystems and living 
organisms [2-4] . One of the most fascinating phenomena in the behavior of complex dynamical 
systems made up of many elements is the spontaneous emergence of order and the phenomenon 
of collective synchronization [5], where a large number of the system's constituents forms a 
common dynamical pattern, despite the intrinsic differences in their individual dynamics. 
Of recent interest are a plenty of biological examples that have become accessible at least 
numerically with the advent of modern computers [6] . 

On the other hand, it has been recently shown that many biological [7,8], social [9], 
and technological [10] systems exhibit an intricate pattern of interconnections in the form of 
complex networks [6]. This structural complexity cannot be described by the couplings of a 
regular network. In order to characterize topologically these complex networks, one computes 
the probability, P{k), that any given element of the network has k connections to other nodes. 
Interestingly, many real- world networks such as the Internet, protein interaction networks and 
social webs [11] can be well approximated by a power-law connectivity distribution, P{k) ~ 
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A;"'*'. Additionally, they are characterized by the existence of key nodes which drastically 
reduce the average distance between nodes, the so-called small- world property [12]. 

Most of th(^ studios performed so far have scrutinized the structure of complex networks 
and studied prototype models ran on top of these networks [13]. Their peculiar topological 
properties have been shown to lead to radical changes when dynamical processes such as 
epidemic spreading and percolation phenomena are studied on top of complex heterogeneous 
networks [14-17]. However, for biological and other applications, it would be relevant to 
consider the nodes of a given network as dynamical systems with own identity. Examples 
of such dynamical systems are ensembles of coupled and pulse-couple oscillators with and 
without time delay, widely used because of their relevance to natural systems such as chirping 
crickets and flashing fireflies, among others [18,19]. Besides, there are several studies where 
the conditions for complete synchronization in complex networks are scrutinized [20-22]. 

In this paper, we numerically study the synchronization of coupled phase oscillators follow- 
ing Kuramoto's model [23, 24] on the underlying topology of an scale-free (SF) network. We 
report on the system dynamics by computing the conventional order parameter. The onset of 
synchronization is found at a nonzero value of the coupling strength. Remarkably, the tran- 
sition from a desynchronized to a synchronous state can be characterized with the mean-field 
exponent found for globally-coupled oscillators. The heterogeneous character of the network 
allows us to explore the robustness of the synchronized state under a single perturbation as a 
function of the connectivity of the perturbed oscillator. Interestingly, we found that the more 
connected a node is, the more stable it is. 

Let us consider an SF network where each node i (i = 1, • • • , TV) is a planar rotor character- 
ized by an angular phase, 9i, and a natural or intrinsic frequency oji. Two oscillators interact 
if they are linked together by an edge of the underlying network. The individual dynamics of 
the zth node is described by 



where J(i) is the set of neighbors of the rotor i as dictated by the architecture of the network 
and A is the coupling strength, identical for all edges. The natural frequencies and the initial 
values of 6i are randomly drawn from a uniform distribution in the interval (—1/2, 1/2) and 
(— 7r,7r), respectively. On the other hand, in order to produce SF networks wc have used 
the BA procedure [25]. In this model, starting from a set of toq nodes, one preferentially 
attaches at each time step a newly introduced node to m older ones (m = 3 has been set). 
The procedure is repeated many times and a network with a power law degree distribution 
P{k) ~ k~'' with 7 = 3 and average connectivity {k) = 2m = 6 builds up. This network 
is a clear example of a highly heterogeneous network because the degree distribution has 
unbounded fluctuations when N ^ oo. 

The original Kuramoto model corresponds to the simplest case of globally coupled (all- 
to-all), equally weighted oscillators where the coupling strength A = K/N to ensure that the 
model is well behaved in the thermodynamic limit [23,24]. The onset of synchronization occurs 
at a critical value of the coupling strength, = 2/Trg{(jJo), where g{uj) is the distribution 
from which the natural frequencies are drawn evaluated at the mean frequency wq- The 
second-order phase transition is characterized by the following order parameter 




(1) 




(2) 



Y. Moreno et al: Synchronization of Kuramoto Oscillators in SF Networks 3 



1.0 



0.8 - 



0.6 



0.4 - 



0.2 



1.0 p 

0.8 - 

0.6 

0.4 - 

0.2 - 

0.0 ^ 
0.0 



e oN=lO^ 

□ N=10 

o o N=5xlO'' 

□ 



0.0 



0.1 



0.2 



0.0 0.2 0.4 0.6 



0.8 



1.0 



Fig. 1 - Coherence r as a function of A for several system sizes. The onset of synchronization occurs 
at a critical value Ac = 0.05(1). Each value of r is the result of at least 10 network realizations and 
1000 iterations for each A*'. The inset is a zoom around Ac. 



which behaves when both N oo and < — > oo as r {K - for K>Kc being /3 = 1/2. 

In order to inspect the dynamics of the N oscillators on top of complex topologies, we have 
performed extensive numerical simulations of the model in BA networks. Starting from A = 0, 
we increase at small intervals its value. Then, we integrate the equations of motion Eq. (|T]l 
over a sufficiently large period of time (at least 10^ integration steps) to ensure that the system 
is in a stationary state, and the order parameter r is computed. The procedure is repeated 
gradually increasing A until the system evolves to a state of collective phase synchronization. 

In the case of random SF networks the global dynamics of the system is qualitatively the 
same as for the original Kuramoto model as shown in Fig. ^ for several system sizes. As the 
coupling is increased from small values, the strength of the interactions is not enough to break 
the incoherence produced by their individual dynamics. This behavior persists until a certain 
critical value Ac is crossed. At this point some elements lock their relative phase and a cluster 
of synchronized nodes comes up. This constitutes the onset of synchronization. Beyond this 
value, the population of oscillators is split into a partially synchronized state made up of 
oscillators locked in phase that adds to r and a group of nodes whose natural frequencies are 
too spread as to be part of the coherent pack. Finally, after further increasing the value of 
A, more and more nodes get entrained around the mean phase and the system settles in a 
completely synchronized state where r « 1. 

This picture is clearly illustrated in Fig. |3 where we have plotted the distributions of 
phases D{9) for five different values of A between and 1. Below the critical point the phases 
are uniformly scattered through the entire interval (— 7r,7r). When A > Ac the distribution 
shrinks around the mean value 9 = and the dispersion gets smaller as A grows, signaling 
that the system is in a synchronous state. 
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Fig. 2 - Normalized phase distributions D{6) for different values of the control parameter A. The 
curves depicted correspond to values of A below, near and above Ac as indicated. The network is 
made up of A'' = 10* nodes. 

Now we proceed to investigate the critical parameters of the system dynamics. First, we 
should determine with precision the exact value of Ac. This is not an easy task because there 
are several sources of heterogeneity and averages should be properly taken —from network 
realizations to the initial distributions of 6i and Ui— through lengthly numerical calculations. 
Additionally, finite-size effects come into play. We have determined Ac and studied the critical 
behavior near the synchronization transition by means of a standard finite-size scaling (FSS) 
analysis [26]. For a given network size N, we have that no synchronization is attained below Ac 
and that r{t) decays to a small residual value of size 0{1/^/N). Hence the critical point may 
be found by examining the A^-dependence of r{\,N). For A < Ac (sub-critical regime), the 
stationary value of r falls off as iV~^/^, while for A > Ac, r approaches a nonzero, stationary 
value as — > oo though still with 0{1/^/N) fluctuations. In this way, plots of r versus N as 
that of Fig. allow us to locate the critical point Ac- 

Following the FSS procedure, our best estimate gave a value for the critical coupling 
strength Ac = 0.05(1). Besides, we found that r ^ (A — Ac)'' above the critical point with 
P = 0.46(2) indicating that the square-root behavior typical of the mean-field version of the 
model (all-to-all architecture) also holds for SF networks. In addition, it is worth stressing 
the very existence of a critical point. This is the opposite to what has been found when 
percolation or epidemic spreading models are ran on top of complex heterogeneous networks 
even if correlations are taken into account [14-17,27]. Furthermore, the critical point shifts 
to the left when the average connectivity (fc) of the underlying network increases, but it 
is always distinctly different from zero. Other numerical calculations (not shown) indicate 
that the relation A^''^ • (k) is roughly constant when (k) varies, which implies that there is a 
non-trivial critical point even in the infinite size limit [28] . 
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Fig. 3 - Finite-size scaling analysis used to determine the value of the critical point. The curves 
suggest the existence of a critical point above 0.04. The values of r are collected after 10* integration 
steps and averaged over a time window of 100 additional integration steps. Finally, each point is 
averaged over at least 10 network realizations and 10'^ different initial conditions. 



Once we have characterized the emergence of spontaneous synchronization, we look at the 
stability and robustness of the synchronous state. The most interesting and influential topo- 
logical property of complex heterogeneous networks is that the fluctuations of the connectivity 
distribution are unbounded. In other words, there is a clear distinction between the nodes ac- 
cording to their connectivities. From a practical point of view —and worth taking into account 
as a design principle in natural or artificial networks— , it would be particularly relevant that 
the most highly connected nodes were also the most robust with respect to synchronization 
when they were perturbed. We have computed the average time (t) it takes for a node to 
be again in the synchronized cluster as a function of its connectivity k after being perturbed 
and put out of synchronization. This is done by assigning to a randomly choosen node i of 
connectivity ki a new phase 9i which differs in — tt to its synchronization value. In this way, 
all nodes are perturbed in the same amount and one may calculate the time it takes for i to 
recover from the perturbation. Note that as A is high enough, the oscillator i ends up in the 
synchronous state with the same 9i it had before being perturbed. The results are drawn in 
Fig. 01 where a clear power-law dependency (r) ~ k~'^ , with v = 0.96(1), can be identified. 
Hence, the more connected a node is, the more stable it is. The power-law behavior points to 
an interesting result, namely, it is more easy for an element with high k to get locked in phase 
with its neighbors than for a node linked to just a few others. Furthermore, the destabilization 
of a hub does not destroy the synchrony of the group it belongs to. Instead, it works the other 
way around, the group formed by the hub's neighbors recruits it again. 

This behavior and the dependency (r) ^ k~^ may be understood by the following simple 
argument. As we are perturbing a single node i and this perturbation, is small, we can 
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Fig. 4 - Log-Log plot of the dimensionless average time (r) it takes for a node of connectivity k to 
be back in the synchronous state after being perturbed. The least square fit to the data gives for the 
exponent v = 0.96(1). The results were averaged over 10 network realizations and 500 perturbations 
for each k. X is set to 0.4. See the text for details on the definition of r. 



assume that it affects only the first neighbors of the perturbed node. Hence, the stability 
analysis can be locally reduced to the problem of how such a perturbation relaxes in a star 
topology (a single perturbed hub attached to fc ^ 1 oscillators). This approximation is 
particularly suitable in random networks with arbitrary degree distributions (like the BA 
one) because the probability of finding loops (tringles, cycles, etc) is small and vanishes as 
N grows. Linearization of Eq. for this configuration leads to = £,i{0)e'^'*, where rji is 
the eigenvalue corresponding to the oscillator i. Henceforth, the times (t)'s are given by the 
inverse of the eigenvalues, which for a star configuration are rji = —1, for i = 1, ■ ■ ■ , N — 2 
and riN-i = rjhub = —N = —khut — 1 [29]. In other words, the fastest relaxation rate in a 
start topology corresponds to the hub and goes fike 1/khub for khub ^ 1, while the rest of the 
oscillators all have the same relaxation times. Additionally, the eigenvector associated to the 
eigenvalue of the hub indicates that it moves a lot while the other nodes change very little. 
Finally, if we superpose the effects of many perturbations to different nodes (each one being a 
khub + 1 star), it comes out that each of them contributes with 1/khub, khub = kmim ■ ■ ■ , kmax 
to the (r) dependency with k, leading to the law (r) ^ k~^ depicted in Fig.^J 

In summary, we have studied the synchronization of Kuramoto oscillators on top of complex 
scale- free networks. We have found that the onset of synchronization occurs at a nonzero value, 
though small, with a critical exponent around 0.5. We also found that when the synchronous 
state has been attained, highly connected nodes are more robust under perturbations and they 
are recovered in a time which depends on their degree as a power law with exponent close to 
— 1. This law support the suggestion that the actual topology of scale- free networks may be 
a result of some kind of optimization mechanism at a local scale, optimizing in this case the 
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fitness for synchronization of highly connected nodes. This question as weh as the influence of 
time delay and noise on the synchronization of complex networks make it necessary the study 
in more details the connection between graph theory, dynamics on networks and nonlinearity 
in future modeling of complex networks. 
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